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Abstract 

Results of a high-statistics, multi- volume Lattice QCD exploration of the deuteron, the di-neutron, 
the H-dibaryon, and the system at a pion mass of m,r ~ 390 MeV are presented. Calculations 

were performed with an anisotropic nj = 2 -|- 1 Clover discretization in four lattice volumes of 
spatial extent L ~ 2.0,2.5,3.0 and 4.0 fm, with a lattice spacing of hg ^ 0.123 fm in the spatial- 
direction, and bt ~ 6^/3.5 in the time-direction. The is found to be bound by 5=-=- = 
14.0(1.4)(6.7) MeV, consistent with expectations based upon phenomenological models and low- 
energy effective field theories constrained by nucleon-nucleon and hyperon-nucleon scattering data 
at the physical light-quark masses. We find weak evidence that both the deuteron and the di- 
neutron are bound at this pion mass, with binding energies of = 11(05)(12) MeV and Bnn = 
7.1(5.2)(7.3) MeV, respectively. With an increased number of measurements and a refined analysis, 
the binding energy of the H-dibaryon is Bh = 13.2(1.8) (4.0) MeV at this pion mass, updating our 
previous result. 



1 



Contents 



I. Introduction I3] 

II. Lattice QCD Calculations H 

A. Liischer's Method for Two-Body Systems Including Bound States H] 

B. Computational Overview [7] 

III. Baryons and Their Dispersion Relations [7] 

IV. Two-Body Bound States M 

A. The Deuteron M 

B. The Di-Neutron [13 

C. The H-Dibaryon M 

D. E-E- m 

E. s-s- m 

V. Conclusions ED 

References 123] 



2 



I. INTRODUCTION 



A major objective for nuclear physicists is to establish the technology with which to reliably 
calculate the properties and interactions of nuclei and to be able to quantify the uncertain- 
ties in such calculations. Achieving this objective will have broad impact, from establishing 
the behavior of matter under extreme conditions such as those that arise in the interior of 
neutron stars, to refining predictions for the array of isotopes produced in nuclear reactors, 
and even to answering anthropic questions about the nature of our universe. While nuclear 
phenomenology generally describes experimentally measured quantities, its ability to make 
high precision and accurate predictions for quantities that cannot be accessed experimen- 
tally is limited. This situation is on the verge of dramatically improving. The underlying 
theory of the strong interactions is known to be quantum chromodynamics (QCD), and 
the computational resources now available are beginning to allow for ab initio calculations 
of basic quantities in nuclear physics. With further increases in computational power and 
advances in algorithms, this trend will continue and our understanding of, and our ability 
to calculate, light and exotic nuclei will be placed on a solid foundation. 

In nature, two nucleons in the ^5*1 —^Di coupled channels bind to form the simplest 
nucleus, the deuteron {J'" = 1+), with a binding energy of = 2.224644(34) MeV, and 
nearly bind into a di-neutron in the ^5*0 channel. However, little is known experimentally 
about possible bound states in more exotic channels, for instance those containing strange 
quarks. The most famous exotic channel that has been postulated to support a bound state 
(the H-dibaryon [T]) has the quantum numbers of AA (total angular momentum J'^ = 0~^, 
isospin 1 = and strangeness s = —2). In this channel, all six quarks in naive quark models, 
like the MIT bag model, can be in the lowest-energy single-particle state. Additionally, more 
extensive analyses using one-boson-exchange (OBE) models |2J and low-energy effective field 
theories (EFT) [Sj H] , both constrained by experimentally measured nucleon-nucleon (NN) 
and hyperon-nucleon (YN) cross-sections and the approximate SU(3) flavor symmetry of 
the strong interactions, suggest that other exotic channels also support bound states. In the 
limit of SU(3) flavor symmetry, the ^S'o-channels are in symmetric irreducible representations 
of8(g)8 = 27©10©T0©8©8©1, and hence the E'E', S'S-, and nn (along with 
nS~ and S^S^) all transform in the 27. YN and NN scattering data along with the leading 
SU(3) breaking effects, arising from the light-meson and baryon masses, suggest that E^E~ 
and E~E~ are bound at the physical values of the light-quark masses |2], |3l H]. 

Recently, the first steps have been taken towards calculating the binding energies of light 
nuclei directly from QCD. Early exploratory quenched calculations of the NN scattering 
lengths |5], |6] performed more than a decade ago have been superseded hj Uf = 2 + 1 
calculations within the last few years [H [8] (and added to by further quenched calculations [Q] 
[TU] ^). Further, the first quenched calculations of the deuteron [12], ^He and ''He 
have been performed, along with Uf = 2-1-1 calculations of ^He [13] and multi-baryon 
systems containing strange quarks [H]. Efforts to explore nuclei and nuclear matter using 
the strong coupling limit of QCD have led to some interesting observations [I5]. Recently, 
Uf = 2 + 1 calculations by us (NPLQCD) [16], and subsequent n/ = 3 calculations by the 
HALQCD collaboration [17j, have provided evidence that the H-dibaryon (with the quantum 



The HALQCD collaboration has produced non-local, energy-dependent, and sink-operator dependent 
quantities from lattice spatial correlation functions that contain the same, but no more, information than 
the NN energy eigenvalues in the lattice volume(s), e.g. Ref. pTj . 
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numbers of AA) is bound at a pion mass of m,r ~ 390 MeV [NPLQCD] and m.^ ~ 837 MeV 
[HALQCD] ^. Extrapolations to the physical light-quark masses suggest a weakly bound 
H-dibaryon or a near threshold resonance exists in this channel [THi Il9j . 

In this work, which is a continuation of our high-statistics Lattice QCD explorations [HI [HI 
[2ni EI], we present evidence for S^H^(^S'o) and H-dibaryon (refining of our results presented 
in Ref. [IE]) bound states, and weak evidence for a bound deuteron and di-neutron at a 
pion mass of ~ 390 MeV. The results were obtained from four ensembles of = 2 -|- 1 
anisotropic clover gauge-field configurations with a spatial lattice spacing of bg ~ 0.123 fm, 
an anisotropy of ~ 3.5 and with cubic volumes of spatial extent L ~ 2.0, 2.5, 3.0 and 
4.0 fm. 

In section |Tl| a concise description of the specific LQCD technology and computational 



details relevant to the present two-body bound state calculations are given. Section III 



presents the results of the LQCD calculations of the single baryon masses and dispersion 
relations (critical for understanding bound systems), and in section IV the results for the 
bound states are presented. Discussions and our conclusions can be found in section IVj 



II. LATTICE QCD CALCULATIONS 

Lattice QCD (LQCD) is a technique in which space-time is discretized into a four- 
dimensional grid and the QCD path integral over the quark and gluon fields at each point 
in the grid is performed in Euclidean space-time using Monte Carlo methods. A LQCD 
calculation of a given quantity will differ from its actual value because of the finite volume 
of the space-time (with x T lattice points) over which the fields exist, and the finite 
separation between space-time points (the lattice-spacing). However, such deviations can be 
systematically removed by performing calculations in multiple volumes with multiple lattice 
spacings, and extrapolating using the theoretically known functional dependences on each. 
In the following subsections, we review the details of LQCD calculations relevant to the 
current work and introduce the ensembles studied herein. 



A. Liischer's Method for Two-Body Systems Including Bound States 

The hadron-hadron scattering amplitude below the inelastic threshold can be determined 
from two-hadron energy levels in the lattice volume using Liischer's method [221 [231 121|- In 
the situation where only a single scattering channel is kinematically allowed, the deviation 
of the energy eigenvalues of the two-hadron system in the lattice volume from the sum of 
the single-hadron energies is related to the scattering phase shift, S{k), at the measured 
two-hadron energies. For energy eigenvalues above kinematic thresholds where multiple 
channels contribute, a coupled-channels analysis is required as a single phase shift does not 
parameterize the S-matrix. Such analyses can be performed, but they are not required in 
the current context. The energy shift for two particles A and B, AE = E^b ~ — Eb, can 
be determined from the correlation functions for systems containing one and two hadrons. 



^ One should note that both calculations were performed at approximately the same spatial lattice spacing 
of 0.12 fm. 
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For baryon-baryon systems, correlation functions of the form 



(1) 



Csi,s.;r(pi,P2,i) = E e^P-^^e^P--^r°i^,^(Bi,«,(xi,i)B2,a.(x2,i)Bi,^,(xo,0)B2A(xo,0)), 



XI, X2 



are used, where B denotes a baryon interpolating operator, ctj and f3i are Dirac indices, and 
the r are Dirac matrices that typically project onto particular parity and angular momentum 
states. The (...) denote averaging over the gauge-field configurations and Xq is the location 
of the source. The interpolating operators are only constrained by the quantum numbers of 
the system of interest, and the simplest forms are 

p«(x,t) = e^^*=<(x,t)[M^"^(x,t)C75ci'(x,t) 

A„(x,t) = e'^X(x,t) [M^'^(x,t)C75ci*^(x,t)' 

E+(x,i) = e^^Vjx,^) [u^^{^,t)C^,s''{^,t) 

S°(x,i) = e^^Vjx,t) [u^^iK,t)C^,s\K,t)] , (2) 

where C is the charge-conjugation matrix and ijk are color indices. Other hadrons in the 
lowest-lying octet can be obtained from the appropriate combinations of quark flavors. The 
brackets in the interpolating operators indicate contraction of spin indices into a spin-0 
"diquark" . Away from the time slice of the source (in this case t — 0), these correlation 
functions behave as 



C^f (P,^) = E ^!4(P) ^(P) e--^'(-)*, 

n 



(3) 
(4) 



where Eq^\o) = rriA and E^^^^Q) are the energy eigenvalues of the two-hadron system 
at rest in the lattice volume. The quantities Z^J')x (Z^^x) are determined by the overlap of 
the source (sink) onto the n}^ energy eigenstate with the quantum numbers of X. At large 
times, the ratio 



^ninB (P' ~P' ^) t^op 



(iJ), 



^o,abIPJ^o,abIPJ e 



(5) 



decays as a single exponential in time with the energy shift, ^e'^^\q). In what follows, 
only the case p = is considered. The energy shift of the n*^ two-hadron state, 

A£;(^^) = £;^^)(0) -mA-rriB^ \[K + ^a + \/kl + m% -mA-rriB , (6) 
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determines a squared momentum, k'^ (which can be either positive or negative). Below 
inelastic thresholds, this is related to the real part of the inverse scattering amplitude via ^ 

knCOt6{kn) = ^s(kl(^)'] , (7) 



27r. 



where 



l<A 



Six) = hm y — 47r A , (8) 

thereby implicitly determining the value of the phase shift at the energy AE!^^^^ (or the 
momentum of each particle in the center of momentum (CoM) frame, kn), S{kn) |22l [23], 
EU |25l |26]. Thus, the function A;cot5 that determines the low-energy elastic-scattering 
cross-section, A{k) oc {k cot 6{k) — ik)~^, is determined at the energy 

In a channel for which one pion exchange (OPE) is allowed by spin and isospin consider- 
ations, the function k cot 6{k) is an analytic function of |kp for |k| < as determined 
by the t-channel cut in the scattering amplitude. In this kinematic regime, k cot 6{k) can be 
expressed in terms of an effective range expansion (ERE) of the form 

kcot6{k) = + ^ ro |kp + ... , (9) 

where a is the scattering length (with the nuclear physics sign convention) and tq is the 
effective range. While the magnitude of the effective range (and higher terms) is set by 
the pion mass, the scattering length is unconstrained. For scattering processes where OPE 
does not contribute, the radius of convergence of the ERE of k cot 6 is set by the lightest 
intermediate state in the t-channel (or by the inelastic threshold). 

In the situation where a channel supports a bound state, the energy of the bound state 
at rest is determined by eq. ([t]). For fc^^^ < 0, and setting = in, eq. ([t]) becomes 

^ Y — g-l'"!''^ = - 
L ^„ |m| L 



^cot5(A:)|,_ + = I E ^^"1-1^^ = j F^'\nL) , (10) 

where 



fW(/s:L) = 6e-"^ + 6^2 e"^"^ + ^ e"^''^ + ... . (11) 



Perturbation theory can be used to solve eq. (10) when the extent of the volume is much 
larger than the size of the bound system, giving 



72 1 



Ko + ^F(°)(/€oL) + C(e-2"«^/L) with = , ' .(12) 



^ ^ ^ •a-2KoAkcotSl 



^ Calculations performed on anisotropic lattices require a modified energy-momentum relation, and, as a 
result, eq. ([6| becomes 

where £,a,b are the anisotropy factors for particle A and particle B, respectively, determined from the 
appropriate energy-momentum dispersion relation. The masses and energy splitting are given in terms of 
temporal lattice units and fc„ is given in spatial lattice units. 
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kq is the solution to 

kcot5ik)\,^^^^ + Ko = , (13) 

which recovers cot = +i, and is the infinite-volume binding momentum of the 

system. This analysis has recently been extended to bound systems that are moving in the 
lattice volume [271 [28]. 

B. Computational Overview 

Anisotropic gauge field configurations have proven useful for the study of hadronic spec- 
troscopy, and as the calculations required for studying multi-hadron systems rely heavily 
on spectroscopy, considerable effort has been put into calculations with clover-improved 
Wilson fermion actions with an anisotropic discretization. In particular, the Uf = 2 + 1 
fiavor anisotropic Clover Wilson action [291 EO] with two steps of stout-link smearing [31] 
of the spatial gauge fields in the fermion action with a smearing weight of p = 0.14 has 
been used |32i I33j . The gauge fields entering the fermion action are not smeared in the time 
direction, thus preserving the ultra-locality of the action in the time direction. Further, a 
tree-level tadpole-improved Symanzik gauge action without a 1 x 2 rectangle in the time 
direction is used. Anisotropy allows for a better extraction of the excited states as well as 
additional confidence that plateaus in the effective mass plots (EMPs) formed from the cor- 
relation functions have been observed, significantly reducing the systematic uncertainties. 
The gauge field generation was performed by the Hadron Spectrum Collaboration (HSC) 
and by us, and these gauge field configurations have been extensively used for excited hadron 
spectrum calculations by HSC [3l|35l [Ml EH |38l [39] . 

The present calculations are performed on four ensembles of gauge configurations with 
X T of 16^ X 128, 20^ x 128, 24^ x 128 and 32^ x 256 lattice sites, with an anisotropy 
of bt = bs/C, with ^ ~ 3.5. The spatial lattice spacing of each ensemble is bg ~ 0.1227 ± 
0.008 fm, giving spatial lattice extents of L ~ 2.0, 2.5, 3.0 and 4.0 fm respectively. The same 
input light-quark mass parameters, btmi = —0.0840 and btrris = —0.0743, are used in the 
production of each ensemble, giving a pion mass of m-j^ ~ 390 MeV. The relevant quantities 
to assign to each ensemble that determine the impact of the finite lattice volume are mT^L 
and rriTrT, which for the four ensembles are rriTrL ~ 3.86,4.82, 5.79 and 7.71 respectively, and 
m^T ~ 8.82, 8.82, 8.82 and 17.64. 

For the four lattice ensembles, multiple light-quark propagators were calculated on each 
configuration. The source location was chosen randomly in order to minimize correlations 
among propagators. On the {16^ x 128, 20^ x 128, 24^ x 128, 32^ x 256} ensembles, an 
average of {224, 364, 178, 174} propagators were calculated on each of {2001, 1195, 2215, 
739} gauge field configurations, to give a total number of ~ {4.5, 4.3, 3.9, 1.3} x 10^ light- 
quark propagators, respectively. 

III. BARYONS AND THEIR DISPERSION RELATIONS 

The single hadron masses calculated in the four different lattice volumes are given in 
Table |Tj Detailed discussions of the fitting methods used in the analysis of the cor- 
relation functions are given in Ref. [3 [TH |20l HQ]. Infinite volume extrapolations 
of the results obtained from the four ensembles were performed in Ref. [21], and are 
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TABLE I: Results from the Lattice QCD calculations in four lattice volumes with a pion mass of 
rriTT ~ 390 MeV, a spatial lattice spacing of bs ~ 0.123 fm, and with an anisotropy factor of ^ ~ 3.5. 
Infinite-volume extrapolations [21] are shown in the right column. The masses are in temporal 
lattice units (t.l.u). 



L3 X T 


16^ X 128 


20^ X 128 


243 X 128 32^ X 256 


Extrapolation 


L (fm) 


~ 2.0 


~2.5 


~3.0 ~4.0 


00 




3.86 


4.82 


5.79 7.71 


00 




8.82 


8.82 


8.82 17.64 


00 


Mn (t.l.u.) 
Ma (t.l.u.) 
Ms (t.l.u.) 
Ms (t.l.u.) 


0.21004(44)(85) 0.20682(34)(45) 0.20463(27)(36) 0.20457(25)(38) 
0.22446(45) (78) 0.22246(27)(38) 0.22074(20)(42) 0.22054(23)(31) 
0.22861(38)(67) 0.22752(32) (43) 0.22791(24)(31) 0.22726(24)(43) 
0.24192(38)(63) 0.24101(27)(38) 0.23975(20)(32) 0.23974(17)(31) 


0.20455(19)(17) 
0.22064(15)(19) 
0.22747(17)(19) 
0.23978(12)(18) 



shown in the right-most column in Table [Tj In physical units, the extrapolated baryon 
masses are Mn = 1151.3(1. 1)(1.0)(7.5) MeV, Ma = 1241.9(0.8)(1.1)(8.1) MeV, Ms = 
1280.3(1.0) (1.1) (8.3) MeV, and Ms = 1349.6(0.7) (1.0) (8.8) MeV [21]. The difference be- 
tween a mass calculated in a finite lattice volume and its infinite-volume extrapolation is 
due to contributions of the form ~ g-m^ri Such deviations must be small compared to 
the two-body binding energies to ensure that the finite volume bindings are due to the T- 
matrix |4Tl 22] and not from finite volume distortions of the forces. It has been shown [TU| [21] 
that the largest two volumes, the 24^ x 128 and 32^ x 256 ensembles, are sufficiently large 
to render the ~ e"™'"^ modifications to Liischer's eigenvalue relation negligible at the level 
of precision we are currently able to achieve. In what follows, we only consider results from 
these ensembles. 

Liischer's method assumes that the single-hadron energy-momentum relation is satisfied 
over the range of energies used in eq. ([T]). In order to verify that the energy-momentum 
relation is satisfied, single hadron correlation functions were formed with well-defined lattice 
spatial momentum p = for |np < 5. Retaining the leading terms in the energy- 
momentum relation, including the lattice anisotropy ^, the energy and mass of the hadron 
(in temporal lattice units (t.l.u)), and the momentum in spatial lattice units (s.l.u) are 
related by 

{b,E^{\n\')y = {hMnf + ^ |np . (14) 

The calculated single hadron energies (squared) are shown in fig. [l] as a function of |np, 
along with the best linear fit. The extracted values of S,h are given in Table [IT[ and are seen 
to be consistent with each other within the uncertainties of the calculation (the value for 
the nucleon is somewhat larger). These values are used to convert the two-hadron energies 
and energy differences from temporal lattice units into spatial lattice units which are then 
used in the Liischer eigenvalue relation. 
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FIG. 1: The squared energy (in (t.l.u.)^) of the single baryon states as a function of = |np, 

related to the squared-momentum, |pp = (x) I'^l^' calculated with the 32'^ x 256 ensemble. The 
blue points are the results of the LQCD calculations with the inner (outer) uncertainties being the 
statistical uncertainties (statistical and systematic uncertainties combined in quadrature). The red 
curves correspond to the best linear-fits. 



TABLE II: The anisotropy parameter, of each hadron from the 32^ x 256 ensemble. The result 
for the TT is included for purposes of comparison. 





N A 


S 




TT 




3.559(27)(08) 3.465(31)(06) 


3.456(35)(07) 


3.4654(55)(14) 


3.466(13)(02) 


IV. 


TWO-BODY BOUND STATES 









Of the baryon-baryon channels that we have explored at this pion mass, the states that have 
an energy lower than two isolated baryons in both the 24^ x 128 and 32^ x 256 ensembles 
and suggest the existence of bound states are the deuteron, the di-neutron, the H-dibaryon, 
and the While a negative energy shift can indicate either a scattering state with an 

attractive interaction or a bound state, Liischer's eigenvalue relation allows us to distinguish 
between the two possibilities. For a bound system in the large- volume limit, the calculated 
value of the energy splitting (or binding momentum) gives rise to —icotS — ?■ +1. We now 
examine each of these channels. 
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A. The Deuteron 



The deuteron is the simplest nucleus, comprised of a neutron and a proton. At the physical 
light-quark masses its binding energy is -B = 2.224644(34) MeV which corresponds to a 
binding momentum of ~ 45.70 MeV (using the isospin averaged nucleon mass of M^v = 
938.92 MeV). As it is a spin-1 system composed of two spin-| nucleons, its wavefunction 
is an admixture of s-wave and d-wave, but at the physical quark masses it is known to 
be predominantly s-wave with only a small admixture of d-wave induced by the tensor 
{L = S = 2) interaction. 

The EMPs associated with the nucleon and the neutron- proton system in the ^5*1 —^Di 
channel are shown in the left panels of fig. [2] and fig. |3] for the two ensembles. The correla- 
tion functions that give rise to these EMPs are linear combinations of correlation functions 
generated using eq. ([T]) but with different smearings of the sink operator (s). The combina- 
tions of correlation functions have been chosen to maximize the extent of the ground-state 
plateaus ^. Extended plateaus are observed in both the one and two nucleon correlation 



0.6 p 
0.5 
^ 0.4- 

3 

W 0.2- 
0.1 



p 

np(3Si) 



"0 10 20 30 40 
t (tl.u.) 



50 



60 



0.03 
0.02 
^ 0.01 
^ 

(N 

:^ -0.01 

-0.02 
-0.03 







10 



20 30 
t (tl.u.) 



40 



FIG. 2: The left panel shows an EMP of the nucleon and of the neutron-proton system in the 
■^^i —^Di coupled channels calculated with the 24^ x 128 ensemble (in t.l.u.). The right panel 
shows the |kp (in (s.l.u.)^) of the neutron-proton system calculated with this ensemble, along with 
the fits. 

functions. The right panels of fig. |2] and fig. |3] show the binding momentum of each parti- 
cle in the CoM obtained by taking ratios of the two-baryon and single-baryon correlation 
functions. The deuteron binding energies in each volume calculated with LQCD are 

^^L=24) ^ 22.3 ± 2.3 ± 5.4 MeV , = 14.9 ± 2.3 ± 5.8 MeV . (15) 



The known finite- volume dependence of loosely bound systems, given in eq. (10), and the 
perturbative relations that follow, allow for an extrapolation of the results in eq. (15) to the 
infinite- volume limit, as shown in fig. |4| giving 

^{L=oo) ^ ii±5±i2MeV (16) 

^ The EMPs result from a matrix-Prony analysis |20l of multiple correlation functions. In determining the 
binding energies, multi-exponential fitting and generalized pencil of function (GPoF) methods |43l I44j are 
used in addition to Matrix-Prony and provides consistent results in each case. 
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FIG. 3: The left panel shows an EMP of the nucleon and of the neutron- proton system in the 
^Si —^Di coupled channels calculated with the 32^ x 256 ensemble (in t.l.u.). The right panel 
shows the |kp (in (s.l.u.)^) of the neutron-proton system calculated with this ensemble, along with 
the fits. 



where the first uncertainty is statistical and the second is systematic, accounting for fitting, 
anisotropy, lattice spacing and the infinite volume extrapolation. Despite having statistically 
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FIG. 4: Results of the Lattice QCD calculations of —i cot 6 versus |kp/m^ in the deuteron channel 
obtained using eq. ([T]), along with the infinite- volume extrapolation using eq. (10). The inner 
uncertainty of the infinite-volume extrapolation is statistical, while the outer corresponds to the 
statistical and systematic uncertainties combined in quadrature. 



significant binding energies in the two lattice volumes, the exponential extrapolation to the 
infinite volume limit produces a deuteron binding energy with significance at ~ la. If the 
uncertainties of both LQCD calculations were reduced by a factor of two, the significance 
of the extrapolated binding energy would increase to ~ 3a if the central values remained 
unchanged. From the curvature of the results of the LQCD calculations in fig. [4], it is clear 
the both of these volumes significantly modify the deuteron at this pion mass. Calcula- 
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tions in somewhat larger volumes, or of moving systems [28], would significantly reduce the 
uncertainty introduced by the volume extrapolation. 

25 

20 



> 15 
S 10 
5 


-5 



Hf = 2+1 
Experiment 



200 400 600 

rrij^ (MeV) 



800 



FIG. 5: The deuteron binding energy as a function of the pion mass. The black circle denotes 
the experimental value. The blue point and uncertainty results from the quenched calculations of 
Ref. [12j, while the red point and uncertainty (the inner is statistical and the outer is statistical 
and systematic combined in quadrature) is our present Uf = 2 + 1 result. 

Our n/ = 2 + 1 result and the recent quenched (n/ = 0) result of Ref. [12] are shown 
in fig. [5} along with the physical deuteron binding energy. Clearly, the large uncertainty 
of our present result does not provide much constraint on the dependence of the deuteron 
binding energy as a function of the light-quark masses, other than to demonstrate that the 
deuteron is likely bound at ~ 390 MeV, qualitatively consistent with the quenched result 
at ~ 800 MeV [12]. 

A number of groups have attempted to determine how the deuteron binding energy 
(and the binding of other nuclei) varies as a function of the light-quark masses using 
EFT [IHl UHl H?! HE] and hadronic models [39]. Such a variation impacts the constraints 
that can be placed on possible time-variations of the fundamental constants of nature from 
the abundance of elements produced in Big Bang Nucleosynthesis (BBN) (see Refs. [50l [51] 
for recent constraints from BBN). With the exception of the analysis of Ref. [H], both of 
the EFT analyses, which use naive dimensional analysis (NDA) to constrain the quark-mass 
dependent dimension-six operators that contribute at next-to-leading order (NLO) in the 
chiral expansion, and the hadronic models of Ref. [IH] , suggest that the deuteron becomes 
less bound as the quarks become heavier near their physical values. The present LQCD 
calculation at a pion mass of m^r ~ 390 MeV is somewhat beyond the range of applicability 
of the EFT analyses and so cannot be directly translated into constraints on the coefficients 
of local operators with confidence. Further, the uncertainty in our calculation is too large to 
be useful in a quantitative way. Nevertheless, our result conflicts with the trend suggested 
in most of the EFT and model analyses, and further studies are necessary to resolve this 
issue. 
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B. The Di-Neutron 



In nature, the di-neutron {nn ^Sq) is very nearly bound. The unnaturally large scatter- 
ing lengths in the ^S'o-channel indicate that a very small increase in the strength of the 
interactions between neutrons would bind them into an electrically neutral nucleus. If the 
binding was deep enough, it would have profound effects on nucleosynthesis. Analyses with 
NNEFT allow for the possibility of both bound and unbound di-neutrons for light-quark 
masses larger than those of nature, while indicating an unbound di-neutron for lighter quark 
masses |l5lll6lll7j. In contrast, a mo del- dependent calculation indicates that the di-neutron 
remains unbound for all light-quark masses ^9] . 

The EMPs associated with the nucleon and the di-neutron system are shown in the left 
panels of fig. [6] and fig. [7} The di-neutron binding energies extracted from the LQCD 
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FIG. 6: The left panel shows an EMP of the neutron and of the neutron- neutron system calculated 
with the 24^ x 128 ensemble (in t.l.u.). The right panel shows the |kp (in (s.l.u.)^) of the neutron- 
neutron system calculated with this ensemble, along with the fits. 
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FIG. 7: The left panel shows an EMP of the neutron and of the neutron-neutron system calculated 
with the 32^ x 256 ensemble (in t.l.u.). The right panel shows the |kp (in (s.l.u.)^) of the neutron- 
neutron system calculated with this ensemble, along with the fits. 
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calculations are 

^(L=24) 



10.4 ± 2.6 ± 3.1 MeV , 51^=^^) = 8.3 ± 2.2 ± 3.3 MeV 



(17) 



The volume extrapolation of the results in eq. (17) is shown in fig. [8| and results in an 
extrapolated di-neutron binding energy of 

S^^=°°) = 7.1 ± 5.2 ± 7.3 MeV (18) 

where the first uncertainty is statistical and the second is systematic. This result is suggestive 
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FIG. 8: The results of the Lattice QCD calculations of —icot5 versus |kp/m^ in the di-neutron 
channel obtained using eq. Q, along with the infinite- volume extrapolation using eq. (10). The 
inner uncertainty of the infinite- volume extrapolation is statistical, while the outer corresponds to 
the statistical and systematic uncertainties combined in quadrature. 

of a bound di-neutron at this pion mass, but at the present level of precision an unbound 
system is also possible. 

Our n/ = 2 + 1 result and the recent quenched {nj = 0) result of Ref. [12,j are shown 
in fig. |9| Clearly, the large uncertainty of our present result does not provide a significant 
constraint on the binding of the di-neutron as a function of the light-quark masses. However, 
the LQCD results suggest that the di-neutron is bound at quark masses greater than those 
of nature. This has implication for future LQCD calculations as there are likely light-quark 
masses for which the di-neutron unbinds, and hence the scattering length becomes infinitely 
large. This implies that, at some point in the future, LQCD may be able to explore strongly 
interacting systems of fermions near the unitary limit. However, if the deuteron remains 
bound at heavier quark masses, as suggested by the current work, it may not be possible 
to tune the light-quark masses (including isospin breaking) to produce infinite scattering 
lengths in the ^Si —^Di and ^Sq channels simultaneously and hence eliminating the possibility 
of the triton having an infinite number of bound states for such a specific choice of light-quark 
masses (unless the deuteron is also unbound for an intermediate range of quark masses) ^ . 



^ Such bound states would be the manifestation of an infrared renormahzation group Umit cycle in QCD, 
as conjectured by Braaten and Hammer |52) . 
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FIG. 9: The di-neutron binding energy as a function of the pion mass. The blue point and 
uncertainty results from the quenched calculation of Ref. [12j, while the red point and uncertainty 
(the inner is statistical and the outer is statistical and systematic combined in quadrature) is our 
present Uf = 2 + 1 result. 



C. The H-Dibaryon 

The prediction of a relatively deeply bound system with the quantum numbers of AA (called 
the H-dibaryon) by Jaffe pj in the late 1970s, based upon a bag-model calculation, started 
a vigorous search for such a system, both experimentally and also with alternate theoretical 
tools. As all six quarks, uuddss, can be in an s-wave and satisfy the Pauli principle, such 
a channel may support a state that is more deeply bound than in channels with different 
flavor quantum numbers. Reviews of experimental constraints on, and phenomenological 
models of, the H-dibaryon can be found in Refs. [SSI EH ESI ES] • While experimental studies 
of doubly-strange (s = —2) hypernuclei restrict the H-dibaryon to be unbound or to have a 
small binding energy, the most recent constraints on the existence of the H-dibaryon come 
from heavy-ion collisions at RHIC |57], effectively eliminating the possibility of a loosely- 
bound H-dibaryon at the physical light-quark masses. However, the analysis that led to 
these constraints was model-dependent, in particular in the production mechanism, and 
may simply not be reliable. Recent experiments at KEK indicate that a near threshold 
resonance may exist in this channel [58] . 

A number of quenched LQCD calculations [591 EOl ED |62l ESI EJ] have previously searched 
for the H-dibaryon, but without success. Recently, we and the HALQCD collaboration have 
reported results that show that the H-dibaryon is bound for a range of light-quark masses 
that are larger than those found in nature [THl [IT]. At present, neither of these calculations 
are extrapolated to the continuum, with both calculations being performed at a spatial 
lattice spacing of bg ~ 0.12 fm. Chiral extrapolations in the light-quark masses of these two 
LQCD calculations were performed in Refs. [T8l [T9] to make first QCD predictions for the 
binding energy of the H-dibaryon at the physical light-quark masses. ^ 



^ These extrapolations are significantly less reliable (rigorous) than the chiral extrapolation of simple 
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In the absence of interactions, the AA-SA^-EE coupled system (all three have the same 
quantum numbers) is expected to exhibit three low-lying eigenstates as the mass-splittings 
between the single-particle states are (from the 32^ x 256 ensemble) 



2 (Ms - Ma) 
Mh + Mtv - 2Ma 



0.01317(13)(19) t.l.u . 
0.003397(61)(65) t.l.u 



(19) 



However, if the interaction generates a bound state, it is unlikely that a second or third 
state will also be bound, and therefore the splitting between the ground state and the two 
additional states will likely be larger than estimates based upon the single-baryon masses. 
The EMPs associated with the A and the system with the quantum numbers of the AA are 
shown in the left panels of fig. [TO] and fig. [TTj The binding energies extracted from the 
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FIG. 10: The left panel shows an EMP of the A and of the lowest state in the AA-HA^-SS system 
calculated with the 24^ x 128 ensemble (in t.l.u.). The right panel shows the |kp (in (s.l.u.)^) of 
the AA-HA^-SS system calculated with this ensemble, along with the fits. 



quantities (such as hadron masses) calculated with LQCD. While for a deeply bound H-dibaryon with a 
radius that is much smaller than the inverse pion mass it is possible to arrive at a chiral EFT construction 
with which to calculate the light-quark mass dependence of H-dibaryon mass in perturbation theory, the 
same construction would not be valid when the radius becomes comparable to or larger than A 
weakly bound state can only be generated nonpcrturbatively, and consequently the quark-mass dependence 
of the binding energy is nontrivial, as is clear from the analyses in the two-nucleon sector, e.g. Refs. |45L 
SSI 1171 [S5] . As a result, the assumption of compactness of the state made in Ref. |TH] is difficult to justify 
over a significant range of predicted binding energies. Further, the simple polynomial extrapolations in 
Ref. [TB] are meant to provide estimates alone and cannot be used to reliably quantify extrapolation 
uncertainties. 
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FIG. 11: The left panel shows an EMP of the A and of the lowest state in the AA-HA^-SS system 
calculated with the 32'^ x 256 ensemble (in t.l.u.). The right panel shows the |kp (in (s.l.u.)^) of 
the AA-HA^-SS system calculated with this ensemble, along with the fits. 



LQCD calculations are 



^(L-24) ^ ^-j,_g2 ± 0.88 ± 0.68 MeV , B 



(L=32) 
H 



14.5 ± 1.3 ±2.4 MeV 



(20) 



which agree within uncertainties with the values given in our earlier paper [16]. The volume 
extrapolation of the results in eq. (20) is shown in fig. 12, and gives an extrapolated H- 
dibaryon binding energy of 



B 



(L=oo) 
H 



13.2 ± 1.8 ±4.0 MeV 



(21) 



where the first uncertainty is statistical and the second is systematic. In Ref. [16], 
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FIG. 12: The results of the LQCD calculations of —i cot 5 versus |kp/m^ in the H-dibaryon channel 
obtained using eq. ([T]), along with the infinite- volume extrapolation using eq. (10). The statistical 
and systematic uncertainties have been combined in quadrature. 



was assigned a volume extrapolation uncertainty of ±1 MeV. In the present analysis, this 
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systematic uncertainty has been reduced to ±0.3 MeV by working to higher orders in the 
volume expansion [28]. The uncertainty in the energy-momentum relation is unchanged, 
and is estimated to be ±0.6 MeV. The updated result in eq. (21) at ~ 390 MeV and the 
result of the nj = 3 calculation at ~ 837 MeV [T7] are shown in fig. 13, Also shown in 
this figure are two naive extrapolations, one that is linear in m^^ and one that is quadratic 

18] . The extrapolations indicate that the LQCD calculations 



in 771,^, as discussed in Ref. 
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FIG. 13: Extrapolations of the LQCD results for the binding of the H-dibaryon. The left panel 
corresponds to an extrapolation that is quadratic in m^r, of the form BHinin) = Bq + dim^. 
The right panel is the same as a left panel except with an extrapolation of the form 13]^(jti-j^^ — 
Bo+diniT^. In each panel, The blue point and uncertainty results from the n/ = 3 LQCD calculation 
of Ref. [T7], while the red point and uncertainty is our present nj = 2 + 1 result. The green dashed 
vertical line corresponds to the physical pion mass. 

are presently not at sufficiently small quark masses to determine if the H-dibaryon is bound 
at the physical light-quark masses. 



D. 



Experimental information on the hyperon-hyperon interactions in the s < —2 sector does 
not exist, presenting a significant handicap to studies of the composition of neutron star 
matter. As an example of the importance of these interactions, Ref. [66] shows that when 
a strongly attractive SS interaction is used in the Tolman-Oppenheimer-Volkoff equation, 
new stable solutions appear, corresponding to compact hyperon stars with masses similar to 
neutron stars but with smaller radii. From the theoretical point of view, the approximate 
flavor SU(3) symmetry of QCD indicates that a bound state in the channel is likely. 

Phenomenological analyses of NN scattering and YN scattering provide a determination of 
the strength of the interaction for two baryons in the 27 irreducible representation of flavor 
SU(3) that also contains the system. The OBE model developed by the Nijmegen 

group, NSC99 [2J ^, which include explicit breaking of flavor SU(3) symmetry by using the 
physical meson and baryon masses, and chiral EFT [67|, predicts a bound state in the 
channel [H H] at the physical pion mass ^. However, only moderate attraction is obtained 



^ The recently developed extended soft-core models do not yet include the s < —2 sectors. 
® The SS('^S'i) and NN(^S'i) states belong to different irreducible representations (10 and 10, respectively) 
and therefore SU(3) flavor symmetry alone is unable to predict whether an analog of the deuteron in the 
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within a constituent quark model 



A small S S interaction was calculated in the 



20^ X 128 ensemble [8] used in this work but may be subject to significant finite volume 
uncertainties. LQCD calculations performed in the fiavor SU(3) limit [69j, in volumes of 
16^ X 32 with a lattice spacing of bg ~ 0.12 fm and at pion masses of 1014 and 835 MeV found 
attractive interactions in the fiavor singlet t-channel responsible for interactions. 

Our present LQCD calculations provide clear evidence for a bound E~E~ state at a pion 
mass of m^r ~ 390 MeV. The EMPs associated with the S and the system are shown 



in the left panels of fig. 14 and fig. 15 
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FIG. 14: The left panel shows an EMP of the H and of the system calculated with the 

24^ X 128 ensemble (in t.l.u.). The right panel shows the |kp (in (s.l.u.)^) of the system 
calculated with this ensemble, along with the fits. 
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FIG. 15: The left panel shows an EMP of the H and of the H~H~ system calculated with the 
32^ X 256 ensemble (in t.l.u.). The right panel shows the |kp (in (s.l.u.)^) of the system 
calculated with this ensemble, along with the fits. 

The S^S^ binding energies extracted from the LQCD calculations are 

B^JrJ^^ = 11.0 ±1.3 ±1.6 MeV , B^JtJ? = 13.0 ± 0.5 ± 3.9 MeV . (22) 



s — —4 sector exists. 
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The volume extrapolation of the results in eq. (22) is shown in fig. 16, and results in an 
extrapolated binding energy of 



B. 



(L=oo) 



14.0 ± 1.4 ± 6.7 MeV 



(23) 



where the first uncertainty is statistical and the second is systematic. This indicates that, at 
the ~ 2a level, the channel supports a bound state. The fact that the binding energy 
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FIG. 16: The results of the Lattice QCD calculations of —icot6 versus |kp/m^ in the 
system obtained using eq. ([T]), along with the infinite- volume extrapolation using eq. (10). The 
inner uncertainty of the infinite- volume extrapolation is statistical, while the outer corresponds to 
the statistical and systematic uncertainties combined in quadrature. 

calculated in the 24^ x 128 ensemble has k cot 6^ indicates that this volume is significantly 
modifying the bound state, and that calculations in larger volumes, or with non-zero 

total momentum, would refine the volume extrapolation. 

This result and the predictions of OBE models and leading order (LO) EFT are shown 
in fig. 17 It is important to note that the uncertainty (and significance) of the LQCD 
result is comparable to that of the OBE models and EFT results. Further, this result 
demonstrates that LQCD is rapidly approaching the situation where it will provide more 
precise constraints on exotic systems than can be achieved in the laboratory. It will be 
interesting to see whether J- PARC [TOj or FAIR [71] can provide constraints on the s = — 3 
and s = —4 systems, as well as on the possible H-dibaryon [72]. The binding energy in 



eq. (23) provides strong motivation to return to OBE models and EFT frameworks and 



determine the expected dependence on the light-quark masses. 



E. S-S- 

As the S~S~ (^5*0) system is in the 27 irreducible representation of flavor SU(3), it 
is also expected to be bound, but by somewhat less than the E~E~ system. While 
the NSC97a-NSC97f models [2] estimate the S~E~ binding, -Bs-s-, to lie in the range 
1.5 MeV< -Bs-s-^ 3.2 MeV, large and negative scattering lengths are found in the S~S~ 
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FIG. 17: The H^H^ binding energy as a function of the pion mass. The black hne denotes 
the predictions of the NSC97a-NSC97f models [2] constrained from nucleon-nucleon and hyperon- 
nucleon scattering data. The orange line denotes the range of predictions by Miller [3], and the 
green line denotes the leading order EFT prediction by Haidenbauer and Meifiner (HM) The red 
point and uncertainty (the inner is statistical and the outer is statistical and systematic combined 
in quadrature) is our present nj = 2 + 1 result. The OBE model and EFT predictions at the 
physical pion mass are displaced horizontally for the purpose of display. 



channel with LO EFT [73] in the absence of Coulomb interactions and isospin breaking 
(these results exhibit non- negligible dependence on the momentum cut-off). On the other 
hand, the constituent quark model of Ref. [6^ finds strong similarities between the behavior 
of the S~E~ and nn interactions, leading to similar values for the phase shifts. Our LQCD 
calculations in this channel are inconclusive. While the ground state in the 24^ x 128 en- 
semble is negatively shifted, the ground state in the 32^ x 256 ensemble is consistent with 
zero, and thus is consistent with both a scattering state and a bound state. 



V. CONCLUSIONS 



We have performed precise Lattice QCD calculations of baryon-baryon systems at a pion 
mass of ~ 390 MeV in four ensembles of anisotropic Clover gauge-field configurations 
with a spatial lattice spacing of hs ~ 0.123 fm, an anisotropy of ^ ~ 3.5 and cubic spatial 
lattice volumes with extent L ~ 2.0, 2.5, 3.0 and 4.0 fm. These calculations have provided 
evidence, with varying levels of significance, for the existence of two-baryon bound states 



from QCD, which are summarized in Table III Our LQCD calculations were performed 



TABLE III: A summary of the two-body binding energies determined in this work. 





Deuteron 


Di-neutron 


H-dibaryon 




Binding Energy (MeV) 


11(05)(12) 


7.1(5.2)(7.3) 


13.2(1.8)(4.0) 


14.0(1.4)(6.7) 
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at one lattice spacing, bs ~ 0.123 fm, but discretization effects are expected to be small 
as they scale as O (fe^) for the Clover action. Consequently, we do not expect them to 
significantly alter our conclusions. A second lattice spacing is required to quantify this 
systematic uncertainty. 

By far the most significant result is that the H-dibaryon is bound at the 3a level at this 
pion mass, improving on results we have already presented in Ref. [I6]. At the ~ 2a level of 
significance, we find that the system is also bound, which is qualitatively consistent 

with an array of hadronic models and EFT analyses of this system at the physical light-quark 
masses. It is interesting to note that the level of precision of the S^S^ binding from LQCD 
is comparable to the level of precision associated with the phenomenological predictions. 
With increasing computational resources directed at these two-baryon systems, the QCD 
prediction will become more precise and eventually become input for phenomenological 
models and will be used to constrain the coefficients appearing in the effective field theories. 

A major goal of Lattice QCD is to postdict the anomalously small binding energy of the 
deuteron. We have presented evidence for a bound deuteron from QCD, however the ~ la 
level of significance is well below "discovery level", and our result should be considered a 
first step toward a definitive calculation. Nevertheless, it is now unambiguously clear that 
a precise determination of the deuteron binding energy can be performed with sufficient 
computational resources. Our result hints that the deuteron is bound, as does the result of 
a previous quenched calculation, at heavy pion masses, in contrast with phenomenological 
analyses and with EFT predictions. We also find suggestions of a bound di-neutron which are 
far from definitive, but are consistent with the quenched result at a heavier pion mass [12]. 
If this remains the case when the calculation is refined, there are light-quark masses between 
~ 140 MeV and m,r ~ 390 MeV for which the scattering length in this channel would 
be infinite and the system would be scale-invariant at low energies. 

Phenomenology based upon flavor SU(3) symmetry indicates that the S^S^ system 
should be more bound than the S~S~ system, which in turn should be more bound than 
the di-neutron (which is nearly bound) at the physical light-quark masses, as these three 
systems are all members of the same 27 irreducible representation of SU(3). Our results are 
consistent with this, but further work is required before definitive conclusions can be drawn. 

The results of the Lattice QCD calculations presented in this paper, which refine and 
broaden our previous work [16], provide clear evidence for bound-states of two baryons 
directly from QCD. With the suggestion of a deuteron and a bound di-neutron at this heavier 
pion mass, there is compelling motivation to invest larger computational resources into 
pursuing Lattice QCD calculations at light-quark masses, and to perform such calculations in 
multiple volumes and with multiple lattice spacings. It is clear that enhanced computational 
resources will enable calculations of the properties and interactions of nuclei from QCD with 
quantifiable and systematically removable uncertainties. 
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